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ABSTRACT 

Continuing our series of papers on the 3-D structure and accurate distances of Planetary Nebu- 
lae (PNe), we present here the results obtained for the planetary nebula NGC 40. Using data from 
different sources and wavelengths, we construct 3-D photoionization models and derive the physical 
quantitities of the ionizing source and nebular gas. The procedure, discussed in detail in the previous 
papers, consists of the use of 3-D photoionization codes constrained by observational data to derive 
the three-dimensional nebular structure, physical and chemical characteristics and ionizing star pa- 
rameters of the objects by simultaneously fitting the integrated line intensities, the density map, the 
temperature map, and the observed morphologies in different emission lines. For this particular case 
we combined hydrodynamical simulations with the photoionization scheme in order to obtain self- 
consistent distributions of density and velocity of the nebular material. Combining the velocity field 
with the emission line cubes we also obtained the synthetic position- velocity plots that are compared 
to the observations. Finally, using theoretical evolutionary tracks of intermediate and low mass stars, 
we derive the mass and age of the central star of NGC 40 as (0.567 ± 0.06)Mo and (5810 ± 600)yrs, 
respectively. The distance obtained from the fitting procedure was (1150 ± 120)pc. 

Subject headings: planetary nebulae: general, individual (NGC 40) - ISM - methods: numerical 



1. INTRODUCTION 

Planetary nebulae are end products of the evolution of 
stars with masses below 8Mq and as such have great 
importance in many fields in astrophysics, from basic 
atomic processes in solar like stars to distant galaxies. 
Although the general pict ure of planet ary nebulae for- 
mation is well understood f|Kwokl 120081 ) many questions 
remain unsolved, such as the mechanism by which the 
stellar ejecta end up forming the many observed mor- 
phologies. 

In the past PNe have been studied with empirical 
methods and one-dimensional photoionization models. 
This scenario has recently been significantly changed 
with modern computational capabilities and software. 
As we have discussed in the previous papers of this se- 
ries (Monteiro, Schwarz, Gruenwald, & Heathcote 2004; 
Montciro. Schwarz, Gruenwald. Gucnthncr. & Heathcote 
20051 : ISchwarz fc Monteiroll2006l) , precise distances are of 
paramount importance to the study of these objects. To 
this extent our work provides precise, self-consistently 
determined distances for objects with comprehensive, in 
most cases spatially resolved, observational constraints. 
These objects can provide valuable calibration to 
pre-existing distance scales as well as self-consistently 
determined physical and chemical quantities. 

One of the major limitations in our previous works is 
related to the three-dimensionality of both density struc- 
ture and velocity field. In the previous papers of this se- 
ries, even though we had well determined density maps, 
the structures had to be defined by a combination of 
parametric surfaces such as spheres and ellipsoids with 
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assumed density gradients to match the observations. 

It is well known that hydrodynamical numerical sim- 
ulations provide the nebular 2-D and/or 3-D density 
distributions and the velocity fields. Therefore, com- 
plex structures as knots, clumps or asymmetries arise 
naturally from physical processes. This is particu- 
larly important in asymmetric planetary nebulae, as 
the determination of physical parameter may be com- 
promised by simplified toy models. Several numeri- 
cal models have been proposed to explain the nebular 
morphologies, su ch as the two-wind interac tion model 
(llcke et all 1 1991) ■ iets (lAkashi fc SokeJ [20081). magnetic 
fields (jGarcia-Segura et al.l l2005 t). and even binarity, 
which accounts for anisotropic distribution of the AGB 
wind, or for the formation of bipolar jets. 

Actual density distributions obtained from these sim- 
ulations are, in general, directly compared to observed 
maps. However, it is well known that the gas density dis- 
tribution - or even the column density projection along 
the line of sight (LOS) - may differ when compared to 
the observed maps. The difference may arise because 
the observations depend on the projected fluxes of cer- 
tain spectral lines, which in turn depends on the three- 
dimensional ionization structure. In this sense, only com- 
bined calculations of both hydrodynamics and photoion- 
ization processes can provide realistic emission maps. 

In this work we present a self-consistent model fit for 
NGC 40 using this new approach, employing both hy- 
drodynamical and photoionization schemes. In the next 
section we present the current knowledge about NGC 40 
and the observational data used in our analysis. In Sec- 
tion 3 we present the numerical scheme and the hydrody- 
namical simulations for this object. The photoionization 
model is described and the results shown in Section 4. 
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Following we discuss the kinematics of the nebular ma- 
terial combining the synthetic emissivities and the simu- 
lated data in Section 5, followed by the discussions and 
conclusions. 

2. NGC 40 

NGC 40 is a well studied object, classified as a low ex- 
citation PN with a WC8 type central star. Several imag- 
ing studies revealed a bright (slightly elliptical) core, a 
large halo a n d filamentary stru c tures (iChu et al. 119871 : 
Bahckl [19871: IBalick et al. 1 119921: IMeaburn et al.l 11996^ . 
Clegg et al.l (|1983| ) presented optical and lUE spectra for 
the bright region of the nebula. Their results indicated 
that the abundances were typical and did not display 
large variations, expected for a WC8 central star. The 
authors argue that the CIV 1549 emission from the neb- 
ular envelope is too strong to have been produced by 
normal thermal processes sugesting that it is a conse- 
quence of processes related to the wind of the central 
star. 

IMeaburn et al.l (|1996[ ) obtained narrow band CCD im- 
ages as well as long slit echelle spectra and showed that 
the outer halo is moving with 31km s~^, with a turbulent 
motion of 7km s~^. The electron temperature for that 
region was calculated to be = (7400 ± 160)K. The in- 
ner halo line profiles are split by about 50 km s~^ which 
is consistent with the splitting seen in the bright core. 
The main outer clumps seen in narrow band images are 
shown to be kinematically associated with barrel shaped 
core, thus not being the result of a true jet. The authors 
also find a d ynamical age of ^ 4000yr from the expansion 
of the core. iKoesterke et al.l (|1998l ) and the dynamical 
ages 3500±500vr fo r NGC 40. 

In'Sabbadin et all (|2000[ ) the authors present detailed 
tomography of NGC 40, obtained with echelle spec- 
troscopy in many slit positions. By assuming a sim- 
ple, direct position-velocity correlation they reconstruct 
the spatial distribution of H'^, 0++ and A^+. The re- 
sults show a complex density structure with densities 
as high as BOOOcm^^. The high de nsities found are in 
agreement with results obtained by iLeal-Ferreira et all 
(|2010)(LGMR10 hereafter) who obtained values around 
2500cm~'^, although the errors are large due to the 
we ak lines used. The e xpansion velocities determined 
bv lSabbadin"eFan (l2000l) are around 25km s ^. The au- 
thors also point out that abundance gradients are likely 
to be present in the nebula, which is also consistent with 
the abundance maps of LGMRIO, where the authors 
present the first spatially resolved spectroscopic mapping 
of NGC 40. 

In the infrared bands, iHonv et al.l (|200lD reported the 
discovery of 21/^m and 30/xm emission features in the 
spectrum obtained from the ISO observatory. Based on 
the presence of the 21/^m feature, they argue that the 
bulk of the dust in the nebula has been produced during 
a carbon-rich phase before the atmospheres of these stars 
became hydrogen poor. The authors also present an es- 
timation of t he dynani i cal ag e of the object of 5000j/r. 
However. iCohen et al.l (|2002| ) ascribed the 21/im feature 
to noise. 

When considering the distance determinations, 
NGC 40 follows more or less the tendency of large inter- 
val that most PNe determ inations have. From the Acker 
catalog (jAcker et al.lll992D we find <d>= 1019±357pc, 



ranging from 620 to 2070 pc. 

2.1. Observational data 

The observational data used to constrain our model 
parameters in the present work was taken from different 
sources of the literature. We use spatially resolved emis- 
sion line maps as our main constraint in the ionization 
structure and spatial density and temperature distribu- 
tion. The spatially resolved data comes from the work 
of LGMRIO described previously. We reffer the reader 
to that work and references therein for greater detailed 
discussion. 

For the absolute Hp flux we use the value of 
iCarrascoet aII(fT98l . The main reason is that the data 
obtained in LGMRIO was calibrated using only one star 
and the sky was not considered to be in photometric 
conditions. In this sense, the fluxes obtained from the 
LGMRIO data set are good for relativ e intensities only . 
We also use the integrated fluxes of iLiu et all ()2004f ) 
(hereafter LLLB04) to help constrain the model inte- 
grated fluxes. Other works mentioned in the introduction 
were also used but with less weight in the flnal compar- 
isons. 

In the infrared band we used data values for NGC 40 
obtained from the 2MASS survey, the IRAS point source 
catalog and also spectrum from the ISO archivetQ. We 
also use fluxes at 9 and 18 /im obtained with the Infrared 
Camera (IRC), and at 65, 90, 140 and 160 /xm using the 
far-Infrared Surveyor (FIS) der ived from the AKARI All- 
Sky S urvey made available in iPhillips fc Marquez-Lugd 
(I20T1I) . 

In order to model the kinematics of NG C 40 we used as 
constraints the line velocities obtained bv lSabbadin et all 
(|2000( ) , which obtained spatially resolved, long-slit echel- 
lograms of the nebula a.tU,B and V filters and for differ- 
ent position angles. This data is used for the construc- 
tion of synthetic tomographic maps of the nebula that 
are directly compared to the observations. 

3. HYDRODYNAMICAL SIMULATIONS OF NGC 40 

In order to obtain more realistic density and veloc- 
ity distributions for the photoionization models, we per- 
formed a number of 2.5D hydrodynamical simulations of 
planetary nebula ejecta. 

We have employed the grid code Godunov-MHD, 
which solves the gas hydrodynamical equations in con- 
servative form. The effects of radiative cooling are 
calculated implicitly after each time step (see Falceta- 
Gongalves, Lazarian & Houde 2010, Falceta-Gongalves 
et al. 2010a, b). The high-order shock-capturing scheme 
is based on an essentially non-oscillatory spatial recon- 
struction and Runge-Kutta time integration. The discon- 
tinuities that arise in simulations of supersonic expanding 
shells are b etter solved using the HLLC Riemann solver 
(|Londrillo fc D el Zanna 2000). 

3.1. Setup for PN ejection 

As initial condition for this p roblem we use th e basic 
model of two interacting winds (jlcke et al.lll992f ). which 
represents the ejection of the tip-AGB stellar envelope 
over a preset density distribution given by the previous 

1 ISO data tag: ADS/IRSA.Atlas#2010/1007/114440.6120 
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Fig. 1. — Density (up) and column density (bottom) maps for 
an inclination of 20° with respect to the plane of sky . 

wind of the red giant star. To account for the giant- 
phase stell ar wind density p rofile we set the ambient gas 
density as (jlcke et. alJll989t ): 

Pamb-^(-J , (1) 

where 

/ l3cos29-f} _ 1 \ 

m = l-a(^ J. (2) 

The parameter a is related to the density ratio at the 
polar and equatorial directions, while /3 the steepness of 
the density profile with the latitude. The initial temper- 
ature of the environment is defined as lOOK. 

As mentioned previously, in the past two decades it 
has become clear that the two-wind interaction model 
is unable to reproduce the morphologies of most PNe. 
Many alternative models have been proposed including 
other dynamical mechanisms for the ejection and shaping 
of the PNe, such as jets (Sahai & Trauger 1998, Lee & 
Sahai 2003, Lee & Sahai 2004, Akashi & Soker 2008, Lee, 
Hsu & Sahai 2009) and magnetic fields (Garcia-Segura 
1997, Garcia-Segura et al. 1999, ). From the theoretical 
point of view, if the magnetic pressure is of order of the 
kinetic and thermal energy densities we should expect 
a larger degree of collimation of the lobes. A similar 
effect would be expected in a heavy jet scenario. In a 
situation with any of these two mechanisms the resulting 
a and /3 parameters would be slightly different. In this 
work, however, we decided to use the simpler two-wind 
interaction model in order to reduce the number of free 
parameters in modeling the dynamical properties of this 
specific object. 

We performed 2.5D simulations for different sets of a 
and (3 parameters in order to obtain the most similar 
projected density distribution to the observed morphol- 
ogy for NGC 40. We run models varying a from 0.0 to 



0.9, in steps of 0.1, and (3 from 0.0 to 6.0, in steps of 
0.5. Due to computational limitations it was impossible 
to run the models in 3-dimensions. All models were cal- 
culated with a fixed grid of 512 x 256 resolution. The cell 
size is defined as 2.2 x lO^^cm and the simulation box is 
1.15 X lO^^cm wide. The simulations were run up to a 
dynamical time of lOOOOyrs. The structure adopted for 
the nebulae is constrained by the observed density map 
and the observed projected morphology. The best model 
for this particular object was obtained for a = 0.3 and 
fi = 5.0, as described below. 

3.2. Hydrodynamical results 

The initial distribution of density is determined by the 
anisotropic wind from the AGB-star. At the center of 
the cube, the ejection of the second and faster wind 
compresses the environment gas resulting in the nebu- 
lar structure. The initial anisotropy of the pre-nebular 
wind and of the post-ejection wind results in a roundish 
structure with two lobes. In Figure 1 (up) we show the 
density distribution obtained a.t t — 5900yr. At this time 
the size of the structure is ^ 61500 AU wide. The nebular 
peak density is ^ 2800cm~'^ at the roundish structure. 
Here, the visible filaments are formed by Ray leigh- Taylor 
instability that arise as the lighter and warmer plasma 
expands onto the cold and dense structure. 

At both opposing sides along the major axis of the 
nebula, the round structure is opened forming two lobes. 
As the internal wind material propagates outwards it 
reaches the inner shock. Since the expanding gas is 
not isotropic, the streamlines cross the shock surface 
obliquely resulting in a deviation of the flow. The mate- 
rial flowing along the nebular shell then converges at the 
apex of the lobe, in the so-called "shock focused inertial 
confinement, resulting in a dense structure, as seen in 
Figure 1, and eventually in a jet (Frank, Balick & Livio 
1996). At the stage of the nebula in this simulation there 
is no jet formed. The lobes present numerical densities 
of ~ lOOOcm"^, except for the two opposing knots - also 
formed by RT-instability - at the edges of the major axis 
of the structure, which are 50% denser. The ejecta aver- 
age velocity is ~ 25km s~^. 

Besides the main nebula seen with densities > 
lOOOcm"^, a diffuse gas with densities ~ SOOcm""^ sur- 
rounds the whole structure. The inner surface is deter- 
mined by the reverse shock of the main nebula with the 
expanding hot wind of the central source. The white 
dwarf wind blows the hot and low density (< 200cm^'^) 
gas which is, in general, observed in X-rays. 

Despite of the fact that the simulation is performed 
in 2.5D, the structure may be considered of cylindrical 
symmetry. Therefore, it is possible to rotate the physi- 
cal variables around the nebular major axis of symmetry 
in order to obtain a three dimensional distribution. We 
must point out here that when this approximation is used 
the small scale structures, such as those generated by the 
RT instabilities mentioned above, will appear in the col- 
umn density maps as circles/rings instead of knots, as 
would be expected from 3D simulations. For the pur- 
poses or the present work this limitation does not influ- 
ence our results. With the cube of density it is possible 
to integrate it and obtain the column density for several 
LOS's. In Figure 1 (bottom), we illustrate the projected 
column density for an inclination angle of 20° for the 
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major axis with respect to the hne of sight. The column 
density map shows similar features as seen in the density 
map, except for the missing two opposing knots located 
at the end of the major axis of the nebula. The reason 
for that is the small size of these clumps. Interestingly, 
two opposing knots at the major axis are visible in the 
observed maps of NGC 40. As we show further in this 
paper, the synthetic maps obtained after computing the 
photoionization reveal the existence of the knots for low 
ionization lines, such as SII and NIL 

4. PHOTOIONIZATION MODELS FOR NGC 40 

In the present work we used the photoionization code 
Mocassin 3D (version 2.02.54) described in full detail in 
lErcolano et all (|200l . 

The procedure to study NGC 40 is the same adopted 
in previous papers of this series, i.e. we provide a density 
distribution of the nebular plasma and run the Monte- 
Carlo simulation using the Mocassin 3D code for the ra- 
diative transfer, assuming a given luminosity (L^) from 
the central source. The whole process is iterated until 
reasonable agreement with observational constraints is 
achieved. The main improvement of this work with re- 
spect to the previous ones is the use of a self-consistent 
distribution of density from the numerical hydrodynam- 
ical simulations. In the following subsections we discuss 
the input parameters used. 

4.1. The ionizing source 

The choice of the ionizing source in photoionization 
models is far from being the trivial choice of a blackbody 
curve of a given temperature. Recent developments in 
the modeling of stellar atmospheres have provided more 
complex stellar spectra that can be used with physical 
and chemical parameters, such as log{g) and abundances 
of the elements which may be important. 

In the field of PNe central stars (CS) the standard 
choice of models is the grid provided by T. RauchQ where 
the user can find NLTE stellar model atmosphere fluxes 
which cover the parameter range of PNe ionizing sources: 
Teff = 50 - 190 kK, log{g) = 5-9 (cm/sec^) and distinct 
abundances. 

We have explored a set of CS models to obtain a best 
match to the observed PN line intensities in NGC 40. 
The different spectra used are shown in Figure 2. The 
photoionization code was run with the density distri- 
bution obtained from the hydrodynamical simulation, 
shown in Figure 1 (transformed in 3-D). The plots differ 
from log{g), T^ff and chemical abundances, which are 
described in Table 1. We found best fitting models for 
each CS star studied giving more emphasis on reproduc- 
ing the total Hp fiux first and then other emission lines. 
The final step consisted of fine tuning the abundances of 
elements that are more sensitive to the Tef f such as oxy- 
gen, in an attempt to reproduce the observed line fluxes. 

As expected, the spectra seem very similar at the re- 
gion that contributes most for the ionization of H, mak- 
ing it hard to select an ionizing source based only on 
the total Hj3 flux. To pin point the best spectrum we 
must analyze other emission lines produced at different 
regions of the nebula, and by ions sensitive to the ioniz- 
ing source temperature. We point out that the Hell 4686 

^ available at |http://astro.uni-tuebingeii.de/~rauch/| 
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Fig. 2. — Central star ionizing spectra studied in this work. 

emission line is one of the best to constrain T^f / , however 
it is too weak in NGC 40 and could not be used. The 
most relevant differences appear in the oxygen lines as 
can be seen in Table 1, where the last column shows the 
observed values obtained from LLLB04 and LGMRIO. 

From the data presented in Table 1 we infer that the 
best ionizing spectrum was the one with [H:He:CNO] 
0:0.33:0.67, which are typical values for a PG 1159 type 
star (hydrogen-deficient post-asymptotic giant branch 
stars believed to be descendents of [WC] type stars). 
For a very deta i led r eview of PG 1159 stars see 
IWerner fc Herwi j (|2006( ). Obviously, due to computa- 
tional limitations, the grid search for the ionizing source 
was not exhaustive. However we explored the most likely 
alternatives given the current knowledge of PNe ionizing 
sources and observation data for NGC 40. 

4.2. Dust 

Apart from the main emission lines, NGC 40 also 
presents large continuum fluxes that are related to the 
dust component, as observed in infrared bands. There- 
fore, we have also included dust in the modeling, taking 
advantage of the full potential of the Mocassin 3-D code. 
As discussed previously, some authors have already dealt 
with the infrared data and dust content of this object. 
The main goal in the present work is to establish reason- 
able intervals for the amount of dust present as well as 
some information on its properties, rather than obtaining 
an excellent flt for the infrared observations. 

To constrain the model dust parameters we have gath- 
ered the spectra obtained by ISO as well as photometry 
in the IRAS catalog. The ISO spectra was obtained from 
the NASA/IPAC Infrared Science Archive. We also use 
fluxes at 9, 18, 6 5, 90, 140 and 160 /xm us i ng the data 
made available in iPhillips &: Marguez-Lugj (|2011[ ) . 

The input parameters for the dust obtained from the 
fltting procedure were, assuming the dust to be perfectly 
mixed with the gas, a mass dust to gas ratio of 0.0015 
with a typical MRN size distribution and composition 
of 50% silicates, 40% grafites and 10% PAH grains. The 
data tables for the dust parameters are those available in 
the Mocassin dust data directory and we refer the reader 



3-D structures and distances of NGC 40 



5 



TABLE 1 

Observed and model line fluxes and model central star parameters for NGC 40. 



Observed" 



H:He:CNO 


1:0:0 


0:1:0 


0:1:0 


0:0.33:0.67 


0:0.33:0.67 


0:0.33:0.67 




TefT 


50kK 


50kK 


50kK 


40kK 


50kK 


50kK 




log(g) 


5 


5 


7 


5 


5 


7 




ni3 flux* 


0.032 


0.034 


0.031 


0.033 


0.037 


0.036 


0.038 


He II4686 


0.157 


0.000 


0.000 


0.000 


0.000 


0.000 


0.000 


[N n]5755 


0.051 


0.041 


0.064 


0.03 


0.04 


0.05 


0.03 


[N ii]6548 


0.774 


0.705 


0.852 


0.98 


0.84 


0.90 


0.85 


[N ii]6584 


2.365 


2.154 


2.602 


2.99 


2.56 


2.75 


2.54 


[O II] 3726 


0.078 


0.066 


0.240 


1.110 


0.881 


2.381 


2.571'' 


[0 II] 3729 


0.052 


0.043 


0.149 


0.717 


0.534 


1.442 


1.904'' 


[0 iii]4363 


0.006 


0.005 


0.007 


0.000 


0.004 


0.004 


0.005 


[0 iii]4959 


0.193 


0.170 


0.189 


0.000 


0.172 


0.196 


0.193 


[0 iii]5008 


0.576 


0.509 


0.563 


0.000 


0.514 


0.585 


0.589 


Line Diagnostics 
















Ne: 
















[S ii]6731/6717 


1.166 


1.206 


1.242 


1.236 


1.275 


1.261 


1.308 


[0 ii]3729/3726 
Te: 

[N ii](6584+6548)/5755 


0.667 


0.650 


0.621 


0.640 


0.606 


0.605 


0.739 


61.826 


69.695 


54.138 


132.8 


82.9 


80.7 


74.4 


[O II] 3726/7320 


19.005 


20.575 


16.052 


33.8 


21.7 


21.3 


36.7 


[0 iii](4959+5007)/4363 


126.95 


152.190 


103.041 




171.6 


195.3 


156.4 



relative fluxes from Leal-Ferreira et al. (2010 
relative fluxes from Liu et al. (2004 ) 
1 X W^'^ergcm-^s~^ 



to lErcolano et all 1)20051 ) for details. 

5. PHOTOIONIZATION MODEL RESULTS 

In this section we present the main results from the 
final fitted model for NGC 40. 

5.1. Synthetic observational maps of NGC 40 

Once the emissivity at each wavelength is obtained, 
from the best fit photoionization model, it is possible to 
construct the synthetic observational maps of NGC 40. 
In order to accomplish this we integrated the emissivity 
cube along a given line of sight. In Figure 3 we show the 
synthetic maps obtained for the Ha, [NII]6584, Hell and 
[OIII]5007 spectral lines assuming an inclination of 20° 
for the nebular major axis with respect to the plane of 
sky. 

These synthetic maps cannot be directly compared to 
the observations due to the different spatial field resolu- 
tion. The numerical simulations presented here present 
finer resolution when compared to the observed maps of 
LGMRIO. In this sense, in order to match their spatial 
resolution we convolved the synthetic maps with a 3 arc- 
sec PSF. In Figure 4 we show the original synthetic map 
for the emission line [Nil] 6584 (top) and the convolved 
map for the same line (bottom) . Despite the coarser spa- 
tial resolution, the knots seen at the ends of the major 
axis of symmetry are still visible and the general mor- 
phology of the nebula is kept. At smaller scales, the 
inhomogeneities created by the Ray leigh- Taylor instabil- 
ities are smoothed. 

In Figure 5 we present two density maps obtained from 
the model. The first map in the upper panel was obtained 



from the projection of the electron density data cube 
weighted by the [5//] 6717 emissivity cube. The lower 
panel shows the density map obtained from the usual 
ratio of the convolved [SI I] emission line maps. The two 
show good agreement despite the different methods used 
in the calculation. 

In Figure 6 we present two electron temperature maps 
obtained from the model. The first map in the upper 
panel was obtained from the projection of the electron 
temperature data cube weighted by the [0///]4363 emis- 
sivity cube. The lower panel shows the electron tem- 
perature map obtained from the usual diagnostic ratio 
[O///] 5007-^4959/ [Om] 4363, again using the convolved 
emission line maps. Again, apart from the difference in 
resolution we see good agreement in the values obtained. 

In Figure 7 we present the observed Ha narrow band 
emission map superimposed by the contours of the con- 
volved synthetic map from the best fit model (taking 
into account the ~15% [Nil] contamination in the nar- 
row band filter). The agreement between both, the size 
and morphology of the true nebular emission, and the 
presented model is clear. 

The total line intensities of the best fit model are given 
in Tabled as well as the fitted abundances and ionizing 
source parameters. To estimate a final relative uncer- 
tainty to be used in the fitted model parameters we cal- 
culated the relativer errors for all lines in Table[2l using 
the average of the observed values when both LLLB04 
and LGMRIO are present. The final relative uncertainty 
can be obtained by taking a weighted average of all rel- 
ative errors, using line intensities as the weights. The 
value obtained is 10% which we adopt as the best es- 
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Fig. 3. — Synthetic images obtained from the projection of the data cubes of emissivities calculated by the photoionization code for 4 of 
the most important emission lines. The nebula has been rotated with an angle of 20° with respect to the plane of sky. 



tiniate of the relative la uncertainty of our final fitted 
parameters. 

The model fitting procedure which uses, among other 
constraints, the model image size fitted to the observed 
one for Ha, as well as the absolute H(3 flux, gives a final 
distance of 1150±120pc for NGC 40. 

In Figure 8 we present the total SED obtained from 
the best fit model (solid line), compared to the ISO 
spectra (gray dashed line), UBVRI-JHK (for the cen- 
tral source only) and IRAS photometric data (circles). 
Here, as mentioned before, we assumed the dust to be 
perfectly mixed with the gas and adopted a mass dust 
to gas ratio of 0.0015, with a typical MRN size distri- 
bution and composition. The dust to gass mass ratio 
is in good agreemen t with the value 0.0013 obtained 
IStasinska fc Szczerbal ()1999D . 



6. POSITION- VELOCITY DIAGRAM 

One of the advantages of combining hydrodynamical 
simulations with photionization calculation is the possi- 
bility of obtaining self-consistent PV-diagrams. The syn- 
thetic PV-diagrams for NGC 40 can be di rectly compared 
to th ose obtained observationally by iSabbadin et al.l 

dJooo). 

From the hydrodynamical simulations we ob- 
tained the velocity field (u) and from the pho- 
toionization code M0CASSIN-3D we obtained 
the emissivity for the main nebular spectral 
lines (e^). The line profiles (/„) are ob tained as 
(jFalceta-Gonpalves. Abraham fc Jatenco-Pereira,2006, ): 



^ (27r(T2)l/2 



exp 



(3) 



where the indices % represent each cell of the cube inter- 
cepted by the LOS, ulos is the projected velocity along 
the LOS, t\ is the emission of a given line A at the i-th 
cell and a is the thermal Doppler broadening. 

The result of Equation 3 is a three dimensional vari- 
able PxPxV, that may be converted into a standard 
PV-diagram once a specific direction (the position an- 
gle) is chosen, representing the observational long-slit. 
Detailed observed PV-di a gram s for NGC 40 were pre- 
sented bv lSabbadin et all ( 2000[ ). In this section we com- 
pare their data with the synthetic line profiles obtained 
from the simulations. Both observational (red circles) 
and synthetic (greyscale) data for [N ii] 6584 are shown 
in Figure 9. The synthetic spectral line profiles were ob- 
tained for different inclination angles with respect to the 
LOS, ranging from 6* = 0° (left) to 40° (right), and for 
two position angles for the slit: top - along the minor 
axis of symmetry, and bottom - along the major axis of 
symmetry of the nebula. Visually it is possible estimate 
a best match for 20°. 

In Figure 10 we present the synthetic PV-diagrams 
obtained from the emission maps of H, [Nil], [SII] and 
[OIII], respectively. The circles represent the observed 
data for [N ii] 6584, as comparison. The plots obtained 
for H, [Nil], [SII] are very similar. The knots are visi- 
ble as the darker spots at the edges of the position angle 
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Fig. 4. — Model [Nil] 6584 images obtained from the best fit 
model. Upper image shows the raw model cube projected with 
orientation as described in the text and lower image shows the 
same projection but convolved with a 3 arcsec PSF to match the 
resolution of LGMRIO maps. 

(PA^ ±30°), with a expanding velocity ^ 10km s~^. 
These clumps are not detected in the [OIII] maps, con- 
firming its low ionization. On the other hand, the [OIII] 
PV-diagram obtained along the minor axis of symmetry 
(top), show emission at positions closer to the central 
source, compared to H, [Nil], [SII]. T hese general re- 
sults c orroborate with the observations of lSabbadin et all 
()2000[) , and shows that both the density distribution and 
the velocity field obtained from the simulations agree 
with the observational data of NGC 40. 

7. DISCUSSION AND CONCLUSIONS 

We have obtained a self consistent three dimensional 
model for the PN NGC 40. Unlike previous works we 
have defined the three dimensional density structure by 
using a 2.5D hydrodynamic simulation as described in 
Sec. 3. With the HD density structure we then pro- 
ceeded to the usual 3D photoionization modeling proce- 
dure performed and described in detail in previous papers 
of this series. 

Though the hydrodynamical simulations do include the 
effects of radiative cooling, we do not treat the effects of 
heating and ionization from the central source. However, 
ionization fronts are known to significantly change the 
propagation of shock fronts under certain conditions (see 
Henney et al. 2005 for review). From Henney et al. 2005, 
the role of photoionization advection is related to the 
advection parameter Aad = 7rj~igj defined as the ratio 
between the fluxes of atoms and photons at the ionization 
front, being M the mach number, ^ ~ 12 for typical 
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Fig. 5. — Electron density maps obtained from the photoioniza- 
tion model. Upper panel shows projected electron density cube, 
weighted by the [S//]6717 emissivity cube. Lower panel shows 
electron density obtained from the usual diagnostic ratio using the 
convolved emission line maps. 
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Fig. 6. — Electron temperature maps obtained from the photoion- 
ization model. Upper panel shows projected electron temperature 
cube, weighted by the [0///]4363 emissivity cube. Lower panel 
shows electron temperature obtained from the usual diagnostic ra- 
tio using the convolved emission line maps. 

central star parameters and sound speed at the ionization 
front, and r ^ 10~^^nz, with z being the characteristic 
Stromgren distance. For NGC40, at the evolutionary 
stage analyzed in this work, z ~ 10^''lO^^ and n ~ 3000, 
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TABLE 2 

Observed and model line fluxes and model central star 
parameters for ngc40. 
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we obtain t - 300 - 3000, which leads to Aad ^ 0.004 - 
0.04M, i.e. the thickness of the ionization front is very 
smaU. In such situation the effects of the ionization front 
in the dynamical evolution of the nebula are negligible, 
and the instabilities are in general quenched (Williams 
2002). 

The density and projected column densities obtained 
from the numerical simulations agree with the observed 
large scale morphology of NGC 40. The column density 
projections however do not show the knots at the major 
axis of simmetry of the nebula as observed at low ion- 
ization energy lines, such as [NII]6584. After using the 
simulated nebula as input for the photoionization code 
MOCASSIN 3-D, the resulting emissivities in this line 
showed good agreement with the observed maps also in 
smaller morphological scales such as the knots. This re- 
veals the importance of combining both hydrodynamical 
and photoionization models in order to study the true 
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Fig. 7. — Image comparing the observed Ha narrow band image 
with the contours of the equivalent image from the fitted model. 




A (/nm) 

Fig. 8.— Model SED (soUd Hne) for NGC 40 compared to 
ISO spectra (gray dashed Hne) and UI3VRI, JHK (both for central 
source only) and IRAS photometric data obtained from SIMBAD. 

n ature of the nebular plasma. 

iMartin et al.l (|2002[ ) reported the detection of emission- 
line structures due to the interaction of the PN with the 
interstellar medium (ISM) attributed to Rayleigh- Taylor 
instability. Indeed, this is observed in the simulations. 

In the X-rays band, NGC 40 presents a faint and diffuse 
emission distri buted within a part ial annulus of about 
40" diameter (jMontez et al.l 120051) . The morphology, 
temperature and luminosity inferred from these obser- 
vations indicate that the emission arises from a hot bub- 
ble generated by shocked quasi-spherical fast wind from 
the central star. The results also show no evidence for 
coUimated jets. In the simulations, as shown in Figure 
1, the reverse shock interacts with the continuous wind 
from the central source creating a low density and high 
temperature medium. This region may be associated to 
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Fig. 9. — Position- Velocity diagrams obtained from the emission maps of [Nil] combined witli the velocity cubes from the hydrodynamical 
simulation, for different inclination angles of the PN axis of symmetry and the LOS (0, 10°, 20°, 30°, 40°, respectively). Top row represents 
the position angle along the minor axis of the PN proje cted at the plan e of th e sky, while the bottom row is for the major axis. Circles 
represent observed values extracted from figs. 3 and 4 of lSabbadin et aLl II2000I) . 
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Fig. 10. — Position- Velocity diagrams obtained from the emission maps of H, [Nil], [SII] and [OIII] combined with the velocity cubes 
from the hydrodynamical simulation, for an inclination angle of 20° . Top row represents the position angle along the minor axis of the PN 
projected at th e plane of the sky, wh ile the bottom row is for the major axis. Circles represent observed values for [Nil] extracted from 
figs. 3 and 4 of lSabbadin et al.l gOOO). 



TABLE 3 

Comparison of literature abundance determinations for NGC 40. 



work 


He/H 


O/H 


N/H 


S/H 


C/H 


Perinotto 


4.73e-2 


7.24e-4 


8.86e-5 


2.80e-6 




Liu 


1.2e-l 


4.9e-4 


8.5e-5 


2.6e-6 


6.9e-4 


Clegg 


> 0.044 


8.4e-4 


2.4e-4 


3.9e-6 


l.Oe-3 


Pottasch 


> 0.046 


5.3e-4 


1.3e-4 


5.6e-6 


1.9e-3 


LGMRIO 


1.2e-l 


4.1e-5 








This work 


0.9e-l 


1.9e-4 


6.8e-5 


6.8e-6 


6.5e-4 



the observations in X-rays. chosen as constraints to the modehng. The usual diag- 

The results show good agreement with all observables nostic maps (density, temperature) as well as all emission 
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line maps obtained by LGMRIO are nicely reproduced. 
We also reproduce with a go od degree of agreement the 
kinematical data obtained by iSabbadin et alj ()2000l ) . 

The total line fluxes as well as the total H[3 flux are re- 
produced within the observational uncertainties for the 
fitted distance of 1150 ± 120. The uncertainty in the 
distance was estimated based on the relative errors of 
observed total line fluxes relative to model values and 
estimated errors of observed and model emission line 
map characteristics done by visual comparison. This re- 
sult is in good agreement with the average result and 
uncertainty obtained from many different literature val- 
ues. The good agreement of the size obtained with our 
structure and distance value is clearly seen in Figure 
5, where we overlay the model contours of a simulated 
iJa-|-15%[NII] image with the respective observed one 
from a Ha narrow band filter. 




Log (Teff) 

Fig. 11.— HR diagram for N G C678 1, NGC3132 
I Monteiro. Morisset. Gruenwald. fc ViegasI |[200C|)), NGC6369 
I Monteir o. Schwarz. Gruenwald. Heathcote' {2004')), MZ 1 
I Monteir o. Schwarz, Gruenwald, Gucnthncr, & Heathcote 
1 20051 ) 1, all PNe that had their central star properties de- 
termined by our method. Also plott ed are the literature 

values for comparison. a) IStanghellini, Corra di, Schwarz! 
11993); b) [S tanghcUini et al.| ll2003h c lBae ssgcn ct al. 
|,199d)d)Pottascli & Bcrnard-Sgilagl i201Q)J_ The evolutionary 
tracks arc from Vassiliadis, & Wood (1994'); they are similar to 
the iBlockeri 11995.) models but take metallicities into account. . 

In Figure 11 we compare the values of effective temper- 
ature and luminosity obtained from our model, T^ff — 
(50 ± 5) kK and L = (1719 dz 170) L©, to evolutionary 
tracks of I Vassiliadis, fc WoodI ()1994[ ) and estimate an age 
of (5810 ± 600) yr, for a core mass of (0.567 ± O.O6)M0 
formed from a progenitor with (1 ± 0.1) Mq. The age 
value is in good agreement with literature values men- 
tioned in Section 1. 

The characteristics of the ionizing source and posi- 
tion in the Tgff x L diagram, as shown in Figure 9, are 
consistent with a [WC] star in transition to a PG1159 
type object. This evolutive connection with sequence 



[WCLl^ [WCE]^([WC]-f - PG11 59)^PG1159 was pro- 
posed bv lKoesterke et all (|19980 . It is interesting to note 
the discrepancy between our model results and previous 
ones from empirical methods from the literature. These 
discrepancies are also present in most other objects stud- 
ied using our method. The lower number of assumptions 
and simplifications as well as the large amount of obser- 
vational constraints used in our work indicate that the 
valu es obtained are likely i nore a ccurate . 

In lMorisset &: GeorgievI (|2009l ) the authors argue that 
the method of determining distances with self-consistent 
modeling could be wrong if the existence of a clumping 
factor is assumed. While this is theoretically true, the 
use of clumping factors are not justified in the methods 
developed in our work where the 3D structure used is 
defined down to the cube cell resolution scale. This is 
especially true in the present work since the fiuctuations 
that appear in the density structure we have obtained are 
solely due to hydrodynamical processes and no artificial 
factor is needed. 

In any case, we f ollow the exercise suggested by 
iMorisset fc GeorgievI ()2009[ ) in their Sec. 5.4 to deter- 
mine the clu mping factor. We h a ve ob tained a rela- 
tion from the iVassiliadis, fc WoodI (|1994l ) model tracks 
in the temperature range of the effective temperature 
obtained by our models for NGC 40. The relation 
we obtained was log{age) = 12.025 — 2.543/o(j'(L/LQ). 
It is important to p o int ou t that it is not clear from 
IMorisset fc GeorgievI ([2009) how exactly their relation 
was obtained. We selected points from the evolution- 
ary tracks of the He-burning models in the temperature 
range of 4.6 < log{Teff ) < 4.8 and fitted a straight line 
to the region close to the luminosity range needed. 

Using the luminosity and dynamical age obtained from 
our model fits we find the intersection to the relation 
above and obtain a value of fc 1. Given that this method 
is very crude and indirect and errors are difficult to es- 
timate it is impossible to obtain a reliable confidence 
interval. If we use the lowest literature value for the dy- 
namical age, 3500 years (Koesterke ct al. 1998), we ob- 
tain fc 1.2. Given that the relative error quoted for the 
age is about 15%, we consider this to be consistent with 
a clumping factor of 1, within the errors. 

In line with all the observational constraints used, al- 
though very crude and indirect, the calculations above 
give some independent evidence that there is no need 
for the artificial clumping factor, indicating also the the 
distance determined for NGC 40 is reliable. 

We have also included dust in our model, even though 
we have not explored the added parameter space to fully 
constrain the possibilities. However, with the assump- 
tion that the dust is mixed with the gas and adopting a 
dust to gas ratio by mass of 0.0015 as well as a typical 
MRN composition, we were able to reproduce the main 
characteristics of the observed infrared features as can be 
seen in Figure 6. A discrepancy is seen in the interval of 
10/im to 20/im with the continuum showing lower values 
than observed. The difference could be due to a denser 
inner region (and thus hotter), such as a disk, which we 
have not included. 

In Table 3 we compare the abundances form other 
works to the values obtained from our model fit. The 
first thing to note is the lack of agreement between em- 
pirical literature determinations as well as reliable con- 
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fidence intervals. This is somewhat troubhng since ah 
the empirical methods are quite similar. The major dif- 
ferences are which portion of the nebula was considered 
and which ionization correction factor was used, which 
indicates that much more work needs to be done in these 
areas if a reliable empirical method is to be developed. 
Even so we can say that the resuhs are roughly correct in 
terms of order of magnitude. It is important to note how- 
ever, that our method does not rely on artificial correc- 
tion factors such as filling factors or ionization correction 
factors, leading to results less dependent on assumptions. 
With all the considerations above we believe that our 
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